# Table D.1
rm(list=ls())

(WD <- getwd())
if (!is.null(WD)) setwd(WD)

indata0      = paste0(WD,"/_individual_data/_spfrawdata/","", collapse = NULL)
indata1      = paste0(WD,"/_average_data/_spfrawdata/","", collapse = NULL)

library(haven)
library(AER)
library(sandwich)
library(lmtest)
library(pracma)
library(dplyr)
library(stargazer)
library(plm)
library(pracma)
library(DataCombine)
library(jtools)
library(plyr)
library(ggplot2)


load(paste(indata0,"us_spf_pgdp_ind.Rda",sep=""))
data$ID         = data$id
data            = data[data$qdate>=1970.00,]
data            = data[data$qdate<=2020.00,]
data = ddply(data,.(qdate),transform,f1_ave = mean(f1,na.rm = TRUE))
data = ddply(data,.(qdate),transform,f2_ave = mean(f2,na.rm = TRUE))

data$f1dev      = data$f1-data$f1_ave 
data$f2dev      = data$f2-data$f2_ave 

plm_gold        = plm(f1dev ~ f2dev, data=data,index=c("ID", "qdate"), model="within")
robust_gold     = sqrt(diag(vcovDC(plm_gold, type = "HC1")))


load(paste(indata1,"us_spf_pgdp_avr.Rda",sep=""))
data_avr        = data
data_avr        = data[data_avr$qdate>=1970.00,]
data_avr        = data[data_avr$qdate<=2020.00,]

ls_out_avr       = lm(fc_err ~ fc_rev, data=data_avr)
robust_out_avr   = sqrt(diag(vcovHC(ls_out_avr, type = "HC1")))

stderror = list(robust_out_avr, robust_gold)

stargazer(ls_out_avr, plm_gold, 
          type = "text",
          se = stderror,
          out="_tables/table_d1.txt")